!      	BLOCK DATA FRONTF 
!		include "common.h" 
!		END 

        program XSTREAM 
        implicit none
        integer :: NI,NJ,NK,NIP1,NJP1,NKP1
	double precision :: X(31,31,20),Y(31,31,20),Z(31,31,20),GPJ(31,31,20) 
	double precision :: GQJ(31,31,20),GRJ(31,31,20) 

        NI=30
        NJ=30
        NK=19

        NIP1=NI+1
        NJP1=NJ+1
        NKP1=NK+1
        
!	DO I=1,NIP1
!        DO J=1,NJP1
!        DO K=1,NKP1
!          X(I,J,K)=0.0D0
!          Y(I,J,K)=0.0D0
!          Z(I,J,K)=0.0D0
!        ENDDO
!       ENDDO
!       ENDDO
        
        call UGRID(X,Y,Z) 
        call GRDGEN(X,Y,Z,GQJ,GPJ,GRJ) 
        

! end of XSTREAM 
        contains
!	SUBROUTINE CONCAT(RESLT,STRl,STR2) 
!	CHARACTER *(*) RESLT,STR1,STR2 
!	INTEGER IPOS,IP2 
 
! SEARCH STR1 FOR FIRST EMPTY SPACE 
 
!		IPOS = 1
!		IP2 = 1
!100 	CONTINUE 
!		RESLT(IP2:IP2)=STR1(IPOS:IPOS)
!		IPOS = IPOS + 1
!		IP2 = IP2 + 1
!		IF((STR1(IPOS:IPOS+1).NE.' ').AND.(IPOS.LE.LEN(STR1))) GOTO 100 

!	CONCATENATE THE SECOND STRING

!		IPOS = 1 
!200 	CONTINUE
!		RESLT(IP2:IP2)=STR2(IPOS:IPOS)
!		IPOS = IPOS + 1
!		IF(IPOS.LE.LEN(STR2)) GOTO 200 

!	PAD WITH SPACES 

!300 	CONTINUE
!		RESLT(IP2:IP2)=' '
!		IP2=IP2 + 1
!		IF(IP2.LE.LEN(RESLT)) GOTO 300
!
!		RETURN 
!		END 


	subroutine SAVEGD(X,Y,Z)
!		include "common.h"

!      This subroutine writes the geometric data into a file for later use
!      The file called fname.grd
	double precision :: X(NIP1,NJP1,NKP1),Y(NIP1,NJP1,NKP1),Z(NIP1,NJP1,NKP1)
!        CHARACTER*50 FNAME

!OPEN THE FILE C 
!	CALL CONCAT(FNAME,FILNAM, '.out') 
	open(unit=7,file="grd.out",form='unformatted')
         

!NOW WRITE THE FILE HEADER 
!
	write(7) NI,NJ,NK 
!
!AND THE GRID POINTS 
!
        write(7) X 
        write(7) Y 
	write (7) Z 
!
        close(7)
!
        return
        end subroutine


	SUBROUTINE GRDGEN(X,Y,Z,GPJ,GQJ,GRJ) 
!		include "common.h" 
	INTEGER WHAT,LOOP,I,J,K 
	double precision :: DS,RTDS,TDS,X(NIP1,NJP1,NkP1),Y(NIP1,NJP1,NKP1)
	double precision :: Z(NIP1,NJP1,NKP1),GQJ(NIP1,NJP1,NKP1),GPJ(NIP1,NJP1,NKP1)
	double precision :: GRJ(NIP1,NJP1,NKP1)
	double precision :: accn

	LOOP = 40 
	accn = 0.4D0 
	DS=0.1D0 
	RTDS = 0.0D0 

!      MAIN LOOP STARTS HERE ... 
MAIN_LOOP : DO 
!        CALL SAVEGD(X,Y,Z)
        CALL SAVEGRD(X ,Y,Z)
!        CALL SAVEMTV(X,Y,Z)
	PRINT*,'**** GRID GENERATOR MENU ****' 
	PRINT*,'1 Change number of iterations (currently ',LOOP,' ) ' 
	PRINT*,'2 Change ds (wall distance -currently ',DS, ') ' 
	PRINT*,'3 Gradually decreasing ds '
	PRINT*,'4 Elliptic generator + boundary clustering' 
	PRINT*,'5 Change acceleration parameter' 
	PRINT*,'6 Quit' 
	READ*,WHAT

        SELECT CASE (WHAT)
        
        CASE (1)
        PRINT*, 'Enter no. of iterations ... '
	READ*,LOOP 

        CASE (2)	
	PRINT*, 'Enter ds ...' 
	READ*, DS
	TDS = 0.01D0
	RTDS = 0.005D0 

        CASE (3)
	PRINT*, 'Enter target DS ..'
	READ*,TDS
	PRINT*, 'Enter fractional decrease in ds aer step ..' 
        read*,RTDS

        CASE (4) 
	CALL ELLIPT (X, Y, Z, LOOP, 1.0D0, GPJ, GQJ, GRJ, DS, accn, RTDS,TDS) 

        CASE (5)
	PRINT*, 'Enter accn. '
	READ*, accn 

	CASE (6)
	EXIT MAIN_LOOP
        END SELECT
        END DO MAIN_LOOP
        RETURN
        END SUBROUTINE GRDGEN

        SUBROUTINE ELLIPT (X,Y,Z,LOOP,BC,GPJ,GQJ,GRJ,DS,accn,RTDS,TDS)
!	include "common. h"
        double precision :: X(NIP1,NJP1,NKP1),Y(NIP1,NJP1,NKP1),Z(NIP1,NJP1,NKP1)
        double precision :: BC,accn,xn,yn,zn,xe,ye,ze,xz,yz,zz,xen,yen,zen 
        double precision :: xez,yez,zez,xnz,ynz,znz
        double precision :: g11,g12,g13,g21,g22,g23,g31,g32,g33
        double precision :: GQJ(NIP1,NJP1,NKP1),GPJ(NIP1,NJP1,NKP1),GRJ(NIP1,NJP1,NKP1) 
	double precision :: dummy1,dummy2,dummy3,gap1,Gp,Gq,Gr,IDS,RTDS,TDS,TZRES 
	double precision :: xAterm,yAterm,zAterm,accnx,accny,accnz,DS,XRES,YRES,ZRES 
	double precision :: JAC,a11,a12,a13,a23,a33,a22,dosum,dd,dx,dy,dz,TXRES,TYRES 
	double precision :: AX(100),AY(100),AZ(100),BX(100),BY(100),BZ(100)
	double precision :: CX(100),CY(100),CZ(100),RX(100),RY(100),RZ(100) 
        double precision :: cpu1,cpu2,cputime 
	integer LOOP,L,I,J,K,ITER 
	character*1 STAT 
	STAT='N' 
	accnx = accn 
	accny = accn 
	accnz = accn 
	TXRES = 1.0D0 
	TYRES = 1.0D0 
	TZRES = 1.0D0 
	IDS=DS 

	CALL CPU_TIME(cpu1) 

	DO L=1,LOOP 
	IF ((ABS(TZRES).LT.0.1D0).AND.(RTDS.NE.0.D0)) THEN
	IF (DS.GT.TDS) THEN 
        DS=DS+RTDS*(TDS-DS) 
        ENDIF 
	ENDIF 

	IF (BC.GT.0.0D0) THEN 
	CALL GRDPQ (X,Y,Z,GPJ,GQJ,GRJ,DS,L) 
	CALL SPRDPQ(GPJ,GQJ,GRJ) 
	ENDIF 

	XRES=0.D0
	YRES=0.D0
	ZRES=0.D0
	TXRES=0.D0
	TYRES=0.D0
	TZRES=0.D0

!      K SWEEP from lower wall to upper wall
        DO I=2,NI
        DO J=2,NJ
        DO K=3,NK-1 
        
        xn = (X(I,J+1,K)-X(I,J-1,K))*0.5D0
        yn = (Y(I,J+1,K)-Y(I,J-1,K))*0.5D0
        zn = (Z(I,J+1,K)-Z(I,J-1,K))*0.5D0
        
        xe = (X(I+1,J,K)-X(I-1,J,K))*0.5D0
        ye = (Y(I+1,J,K)-Y(I-1,J,K))*0.5D0
        ze = (Z(I+1,J,K)-Z(I-1,J,K))*0.5D0
        
        xz = (X(I,J,K+1)-X(I,J,K-1))*0.5D0
        yz = (Y(I,J,K+1)-Y(I,J,K-1))*0.5D0
        zz = (Z(I,J,K+1)-Z(I,J,K-1))*0.5D0
        
        xen=(X(I+1,J+1,K)-X(I-1,J+1,K)+X(I-1,J-1,K)-X(I+1,J-1,K))*0.25D0
        yen=(Y(I+1,J+1,K)-Y(I-1,J+1,K)+Y(I-1,J-1,K)-Y(I+1,J-1,K))*0.25D0
        zen=(Z(I+1,J+1,K)-Z(I-1,J+1,K)+Z(I-1,J-1,K)-Z(I+1,J-1,K))*0.25D0
        
        xez=(X(I+1,J,K+1)-X(I-1,J,K+1)+X(I-1,J,K-1)-X(I+1,J,K-1))*0.25D0
        yez=(Y(I+1,J,K+1)-Y(I-1,J,K+1)+Y(I-1,J,K-1)-Y(I+1,J,K-1))*0.25D0
        zez=(Z(I+1,J,K+1)-Z(I-1,J,K+1)+Z(I-1,J,K-1)-Z(I+1,J,K-1))*0.25D0
        
        xnz=(X(I,J+1,K+1)-X(I,J-1,K+1)+X(I,J-1,K-1)-X(I,J+1,K-1))*0.25D0
        ynz=(Y(I,J+1,K+1)-Y(I,J-1,K+1)+Y(I,J-1,K-1)-Y(I,J+1,K-1))*0.25D0
        znz=(Z(I,J+1,K+1)-Z(I,J-1,K+1)+Z(I,J-1,K-1)-Z(I,J+1,K-1))*0.25D0
        
        g11 = yn*zz-yz*zn
        g12 = yz*ze-ye*zz
        g13 = ye*zn-yn*ze
        
        g21 = xn*zn-xz*yn
        g23 = xn*ze-xe*zn
        g22 = xe*zz-xz*ze
        
        g31 = xn*yz-xz*yn
        g32 = xz*ye-xe*yz
        g33 = xe*yn-xn*ye
        
        a11 = g11*g11 + g21*g21 + g31*g31
        a22 = g12*g12 + g22*g22 + g32*g32
        a33 = g13*g13 + g23*g23 + g33*g33
        
        a12 = g11*g12 + g21*g22 + g31*g32
        a13 = g11*g13 + g21*g23 + g31*g33
        a23 = g12*g13 + g22*g23 + g32*g33
        
        JAC = xe*(g11)-xn*(-g12)+xz*(g13)
        
        Gp=GPJ(I,J,K)
        Gq=GQJ(I,J,K)
        Gr=GRJ(I,J,K)
        
        xAterm = -(JAC*JAC)*(Gp*xe+Gq*xn+Gr*xz)
        yAterm = -(JAC*JAC)*(Gp*ye+Gq*yn+Gr*yz)
        zAterm = -(JAC*JAC)*(Gp*ze+Gq*zn+Gr*zz)
        
        !--for x 
        dummy1=-a11*(X(I+1,J,K)+X(I-1,J,K)) 
        dummy2=-a22*(X(I,J+1,K)+X(I,J-1,K)) 
        gap1=-2.0D0*(a12*xen+a13*xez+a23*xnz)
        
        
        !---Assemble TDMA Coefficients for x
        AX(K)=a33
        BX(K)=-2.0D0*(a11+a22+a33)
        CX(K)=a33
        RX(K)=dummy1+dummy2+gap1+xAterm
        
        !---for y 
        dummy1=-a11*(Y(I+1,J,K)+Y(I-1,J,K))
        dummy2=-a22*(Y(I,J+1,K)+Y(I,J-1,K))
        gap1=-2.0D0*(a12*yen+a13*yez+a23*ynz)
        
        !---Assemble TDMA Coefficients for y
        AY(K)=AX(K)
        BY(K)=BX(K)
        CY(K)=CX(K)
        RY(K)=dummy1+dummy2+gap1+yAterm
        
        !---for z 
        dummy1=-a11*(Z(I+1,J,K)+Z(I-1,J,K)) 
        dummy2=-a22*(Z(I,J+1,K)+Z(I,J-1,K)) 
        gap1=-2.0D0*(a12*zen+a13*zez+a23*znz) 
        
        !---Assemble TDMA Coefficients for z
        
        AZ(K)=AX(K)
        BZ(K)=BX(K)
        CZ(K)=CX(K)
        RZ(K)=dummy1+dummy2+gap1+zAterm 
        
        ENDDO 
        !--- end of J loop 
        
        !Implicit BCs for X
        !lower wall
        
        RX(3)= RX(3)-AX(3)*X(I,J,2)
        AX(3)= 0.0D0
        
        !upper wall
        RX(NK-1) = RX(NK-1)-CX(NK-1)*X(I,J,NK)
        CX(NK-1) = 0.0D0 
        
        !Implicit BCs for Y 
        !      lower wall 
        RY(3) = RY(3) - AY(3)*Y(I,J,2)
        AY(3) = 0.0D0
        !upper wall
        
        RY(NK-1) = RY(NK-1)-CY(NK-1)*Y(I,J,NK)
        CY(NK-1) = 0.0D0
        
        !      --Implicit BCs for Z
        !      --lower wall
        RZ(3) = RZ(3)-AZ(3)*Z(I,J,2)
        AZ(3) = 0.0D0
!        upper wall
        RZ(NK-1) = RZ(NK-1)-CZ(NK-1)*Z(I,J,NK)
        CZ(NK-1) = 0.0D0 
        
        !--SOLVE FOR X 
        CALL TDMA(AX,BX,CX,RX,3,NK-1,100) 
        
        !--SOLVE FOR Y 
        CALL TDMA(AY,BY,CY,RY,3,NK-1,100) 
        
        !--SOLVE FOR Z 
        CALL TDMA(AZ,BZ,CZ,RZ,3,NK-1,100) 
        
        
        !      ---Update x,y and z 
        
        DO K=3,NK-1
        
        Z(I,J,K)= accnz*Z(I,J,K)+(1.0D0-accnz)*RZ(K)
        if ((J.eq.2).or.(J.eq.NJ).and.(I.eq.2).or.(I.eq.NI)) then 
        !-- do nothing 
        else if ((J.ne.NJ).and.(J.ne.2)) then
        Y(I,J,K)=accny*Y(I,J,K)+(1.0D0-accny)*RY(K)
        endif 
        
        if((J.eq.2).or.(J.eq.NJ).and.(I.eq.2).or.(I.eq.NI)) then 
        !-- do nothing
        else if ((I.ne.2).and.(I.ne.NI)) then
        X(I,J,K)=accnx*X(I,J,K)+(1.0D0-accnx)*RX(K)
        endif
        
        ENDDO 
        ENDDO 
        ENDDO 
        
        PRINT*,' ITER =',L 
        CALL GRDEXT(X,Y,Z) 
        IF(L.eq.LOOP) THEN 
        
        !--Residual Calculation
        DO K=3,NK-1
        DO J=3,NJ-1
        DO I=3,NI-1 
        
        xn = (X(I,J+1,K)-X(I,J-1,K))*0.5D0
        yn = (Y(I,J+1,K)-Y(I,J-1,K))*0.5D0
        zn = (Z(I,J+1,K)-Z(I,J-1,K))*0.5D0
        
        xe = (X(I+1,J,K)-X(I-1,J,K))*0.5D0
        ye = (Y(I+1,J,K)-Y(I-1,J,K))*0.5D0
        ze = (Z(I+1,J,K)-Z(I-1,J,K))*0.5D0
        
        xz = (X(I,J,K+1)-X(I,J,K-1))*0.5D0
        yz = (Y(I,J,K+1)-Y(I,J,K-1))*0.5D0
        zz = (Z(I,J,K+1)-Z(I,J,K-1))*0.5D0
        
        xen=(X(I+1,J+1,K)-X(I-1,J+1,K)+X(I-1,J-1,K)-X(I+1,J-1,K))*0.25D0
        yen=(Y(I+1,J+1,K)-Y(I-1,J+1,K)+Y(I-1,J-1,K)-Y(I+1,J-1,K))*0.25D0
        zen=(Z(I+1,J+1,K)-Z(I-1,J+1,K)+Z(I-1,J-1,K)-Z(I+1,J-1,K))*0.25D0
        
        xez=(X(I+1,J,K+1)-X(I-1,J,K+1)+X(I-1,J,K-1)-X(I+1,J,K-1))*0.25D0
        yez=(Y(I+1,J,K+1)-Y(I-1,J,K+1)+Y(I-1,J,K-1)-Y(I+1,J,K-1))*0.25D0
        zez=(Z(I+1,J,K+1)-Z(I-1,J,K+1)+Z(I-1,J,K-1)-Z(I+1,J,K-1))*0.25D0
        
        xnz=(X(I,J+1,K+1)-X(I,J-1,K+1)+X(I,J-1,K-1)-X(I,J+1,K-1))*0.25D0
        ynz=(Y(I,J+1,K+1)-Y(I,J-1,K+1)+Y(I,J-1,K-1)-Y(I,J+1,K-1))*0.25D0
        znz=(Z(I,J+1,K+1)-Z(I,J-1,K+1)+Z(I,J-1,K-1)-Z(I,J+1,K-1))*0.25D0
        
        g11 = yn*zz-yz*zn
        g12 = yz*ze-ye*zz
        g13 = ye*zn-yn*ze
        
        g21 = xn*zn-xn*zz
        g23 = xn*ze-xe*zn
        g22 = xe*zz-xz*ze
        
        g31 = xn*yz-xz*yn
        g32 = xz*ye-xe*yz
        g33 = xe*yn-xn*ye
        
        a11 = g11*g11 + g21*g21 + g31*g31
        a22 = g12*g12 + g22*g22 + g32*g32
        a33 = g13*g13 + g23*g23 + g33*g33
        
        a12 = g11*g12 + g21*g22 + g31*g32
        a13 = g11*g13 + g21*g23 + g31*g33
        a23 = g12*g13 + g22*g23 + g32*g33
        
        JAC = xe*(g11)-xn*(-g12)+xz*(g13)
        
        Gp=GPJ(I,J,K)
        Gq=GQJ(I,J,K)
        Gr=GRJ(I,J,K)
        
        xAterm = (JAC*JAC)*(Gp*xe+Gq*xn+GR*xz)
        yAterm = (JAC*JAC)*(Gp*ye+Gq*yn+GR*yz)
        zAterm = (JAC*JAC)*(Gp*ze+Gq*zn+GR*zz)
        
        !-- residual for x 
        dummy1=(X(I+1,J,K)-2.0D0*X(I,J,K)+X(I-1,J,K))*a11 
        dummy2=(X(I,J+1,K)-2.0D0*X(I,J,K)+X(I,J-1,K))*a22 
        dummy3=(X(I,J,K+1)-2.0D0*X(I,J,K)+X(I,J,K-1))*a33
        
        gap1 = 2.0D0*(a12*xen+a13*xez+a23*xnz)
        dosum=(dummy1+dummy2+dummy3)+(gap1+xAterm)
        XRES = XRES + ABS(dosum)
        
        !-- residual for y
        dummy1=(Y(I+1,J,K)-2.0D0*Y(I,J,K)+Y(I-1,J,K))*a11
        dummy2=(Y(I,J+1,K)-2.0D0*Y(I,J,K)+Y(I,J-1,K))*a22
        dummy3=(Y(I,J,K+1)-2.0D0*Y(I,J,K)+Y(I,J,K-1))*a33
        
        gap1 = 2.0D0*(a12*yen+a13*yez+a23*ynz)
        dosum=(dummy1+dummy2+dummy3)+(gap1+yAterm)
        YRES = YRES + ABS (dosum)
        
        !-- residual for z
        dummy1=(Z(I+1,J,K)-2.0D0*Z(I,J,K)+Z(I-1,J,K))*a11
        dummy2=(Z(I,J+1,K)-2.0D0*Z(I,J,K)+Z(I,J-1,K))*a22
        dummy3=(Z(I,J,K+1)-2.0D0*Z(I,J,K)+Z(I,J,K-1))*a33
        gap1 = 2.0D0*(a12*zen+a13*zez+a23*znz)
        dosum=(dummy1+dummy2+dummy3)+(gap1+zAterm)
        ZRES = ZRES + ABS (dosum) 
        ENDDO 
        ENDDO 
        ENDDO 
        
        TXRES =XRES/((NI-1)-2)/((NJ-1)-2)/((NK-1)-2)
        TYRES =YRES/((NI-1)-2)/((NJ-1)-2)/((NK-1)-2)
        TZRES =ZRES/((NI-1)-2)/((NJ-1)-2)/((NK-1)-2) 
        
        if((ABS(TXRES).lt.0.01D0).and.(ABS(TYRES).lt.0.01D0).and.(ABS(TZRES).lt.0.01D0)) then 
                if (STAT.ne.'C') then
                 STAT='C'
                 ITER=L 
        endif
        else
                STAT='N' 
        endif  
                PRINT*, TXRES, ' ' ,TYRES,' ',TZRES
        endif
        PRINT* 
        
        
        dx=X(24,24,3)-X(24,24,2) 
        dy=Y(24,24,3)-Y(24,24,2) 
        dz=Z(24,24,3)-Z(24,24,2) 
        dd=DSQRT(dx*dx+dy*dy+dz*dz) 
        ENDDO 
        
        !-- End of Iteration Loop.
        CALL CPU_TIME(cpu2)
        cputime=cpu2-cpu1
        PRINT*, 'CPU time =',cputime
        PRINT*, 'ITER =' ,ITER 
        !-- check spacing ... 
        
        CALL SEESPC (X, Y, Z) 
        
        RETURN 
        END SUBROUTINE
        
        !-- End of Ellipt 
        
        
        
        
        subroutine TDMA(AA,BB,CC,RR,BL,NL,MAXL) 
        integer :: L,NL,BL,MAXL 
        double precision :: AA(MAXL),BB(MAXL),CC(MAXL),RR(MAXL),TEMP 
        
        !   Forward elimination 
        
        DO L=BL+1,NL 
        TEMP=AA(L)/BB(L-1) 
        BB(L)=BB(L)-CC(L-1)*TEMP 
        RR(L)=RR(L)-RR(L-1)*TEMP 
        ENDDO
        
        !Back substitution
        RR(NL)= RR(NL)/BB(NL) 
        DO L=NL-1,BL,-1 
        RR(L)=(RR(L)-CC(L)*RR(L+1))/BB(L) 
        ENDDO 
        RETURN 
        end subroutine TDMA
        
        
        subroutine GRDPQ(X,Y,Z,GPJ,GQJ,GRJ,DS,L) 
!        inc1ude "common.h" 
        double precision :: X(NIP1,NJP1,NKP1),Y(NIP1,NJP1,NKP1),GQJ(NIP1,NJP1,NKP1),DS
        double precision :: GPJ(NIP1,NJP1,NKP1),Z(NIP1,NIP1,NKP1),GRJ(NIP1,NJP1,NKP1)
        INTEGER I,J,K,L 
        
        K=2
        DO J=2,NJ
        DO I=2,NI 
        CALL SPQR(GPJ,GQJ,GRJ,X,Y,Z,I,J,K,1,L,DS) 
        ENDDO 
        ENDDO 
        RETURN 
        END SUBROUTINE
        
        subroutine SPQR(GPJ,GQJ,GRJ,X,Y,Z,I,J,K,PJ,L,DS) 
!        include "common.h" 
        INTEGER I,J,K,PJ,L 
        double precision :: GPJ(NIP1,NJP1,NKP1),GQJ(NIP1,NJP1,NKP1),X(NIP1,NJP1,NKP1),Y(NIP1,NIP1,NKP1)
        double precision :: Z(NIP1,NJP1,NKP1),GRJ(NIP1,NJP1,NKP1)
        double precision :: xe,ye,ze,xn,yn,zn,xz,yz,zz,DS,sgn,JAC,tmp
        double precision :: xee,yee,zee,xnn,ynn,znn,xzz,yzz,zzz
        double precision :: t1,t2,t3,xen,yen,zen,xez,yez,zez,xnz,ynz,znz
        double precision :: rxe,rye,rze,rxn,ryn,rzn,lxe,lye,lze,lxn,lyn,lzn
        double precision :: rxz,ryz,rzz,lxz,lyz,lzz,xRHS,yRHS,zRHS
        double precision :: g11,g12,g13,g21,g22,g23,g31,g32,g33
        double precision :: rg13,rg23,rg33,lg13,lg23,lg33
        double precision :: a11,a22,a33,a12,a13,a23 
        
        sgn = DBLE(PJ)
        
        xn = (X(I,J+1,K)-X(I,J-1,K))*0.5D0
        yn = (Y(I,J+1,K)-Y(I,J-1,K))*0.5D0
        zn = (Z(I,J+1,K)-Z(I,J-1,K))*0.5D0
        xe = (X(I+1,J,K)-X(I-1,J,K))*0.5D0
        ye = (Y(I+1,J,K)-Y(I-1,J,K))*0.5D0
        ze = (Z(I+1,J,K)-Z(I-1,J,K))*0.5D0
        g13 = ye*zn-yn*ze
        g33 = xe*yn-xn*ye
        g23 = xn*ze-xe*zn
        
        tmp = 1.0D0/g33
        
        zz = DS/DSQRT((g13*tmp)*(g13*tmp)+(g23*tmp)*(g23*tmp)+1.0D0)
        xz = g13*zz*tmp
        yz = g23*zz*tmp
        
        g11 = yn*zz-yz*zn 
        g12 = yz*ze-ye*zz 
        g21 = xz*zn-xn*zz 
        g22 = xe*zz-xz*ze 
        g31 = xn*yz-xz*yn 
        g32 = xz*ye-xe*yz 
        
        a11 = g11*g11 + g21*g21 + g31*g31
        a22 = g12*g12 + g22*g22 + g32*g32
        a33 = g13*g13 + g23*g23 + g33*g33
        a12 = g11*g12 + g21*g22 + g31*g32
        a13 = g11*g13 + g21*g23 + g31*g33 
        a23 = g12*g13 + g22*g23 + g32*g33
        
        JAC = 1.0D0/(xe*(g11)-xn*(-g12)+xz*(g13))
        
        !- second derivatives at lower wall
        tmp = 8.0D0*X(I,J,K+PJ)-7.0D0*X(I,J,K)-X(I,J,K+2*PJ)
        xzz = tmp*0.5D0-sgn*3.0D0*xz
        
        tmp = 8.0D0*Y(I,J,K+PJ)-7.0D0*Y(I,J,K)-Y(I,J,K+2*PJ)
        yzz = tmp*0.5D0-sgn*3.0D0*yz
        
        tmp = 8.0D0*Z(I,J,K+PJ)-7.0D0*Z(I,J,K)-Z(I,J,K+2*PJ)
        zzz = tmp*0.5D0-sgn*3.0D0*zz
        
        xee=X(I+1,J,K)+X(I-1,J,K)-2.0D0*X(I,J,K) 
        yee=Y(I+1,J,K)+Y(I-1,J,K)-2.0D0*Y(I,J,K) 
        zee=Z(I+1,J,K)+Z(I-1,J,K)-2.0D0*Z(I,J,K)
        
        xnn=X(I,J+1,K)+X(I,J-1,K)-2.0D0*X(I,J,K) 
        ynn=Y(I,J+1,K)+Y(I,J-1,K)-2.0D0*Y(I,J,K) 
        znn=Z(I,J+1,K)+Z(I,J-1,K)-2.0D0*Z(I,J,K)
        
        
        !      - xez cross derivatives 
        
        rxe=X(I+1,J,K)-X(I,J,K) 
        rye=Y(I+1,J,K)-Y(I,J,K) 
        rze=Z(I+1,J,K)-Z(I,J,K) 
        rxn=((X(I+1,J+1,K)-X(I+1,J-1,K))+(X(I,J+1,K)-X(I,J-1,K)))*0.25D0 
        ryn=((Y(I+1,J+1,K)-Y(I+1,J-1,K))+(Y(I,J+1,K)-Y(I,J-1,K)))*0.25D0
        rzn=((Z(I+1,J+1,K)-Z(I+1,J-1,K))+(Z(I,J+1,K)-Z(I,J-1,K)))*0.25D0 
        
        rg13 = rye*rzn-ryn*rze
        rg23 = rxn*rze-rxe*rzn
        rg33 = rxe*ryn-rxn*rye
        
        rzz=DS/DSQRT((rg13/rg33)+(rg23/rg33)+1.0D0)
        rxz=(rg13/rg33)*rzz
        ryz=(rg23/rg33)*rzz
        
        lxe=X(I,J,K)-X(I-1,J,K)
        lye=Y(I,J,K)-Y(I-1,J,K)
        lze=Z(I,J,K)-Z(I-1,J,K)
        lxn=((X(I,J+1,K)-X(I,J-1,K))+(X(I-1,J+1,K)-X(I-1,J-1,K)))*0.25D0 
        lyn=((Y(I,J+1,K)-Y(I,J-1,K))+(Y(I-1,J+1,K)-Y(I-1,J-1,K)))*0.25D0
        lzn=((Z(I,J+1,K)-Z(I,J-1,K))+(Z(I-1,J+1,K)-Z(I-1,J-1,K)))*0.25D0 
        
        lg13 = lye*lzn-lyn*lze
        lg23 = lxn*lze-lxe*lzn
        lg33 = lxe*lyn-lxn*lye
        
        lzz=DS/DSQRT((lg13/lg33)+(lg23/lg33)+1.0D0)
        lxz=(lg13/lg33)*lzz
        lyz=(lg23/lg33)*lzz
        
        xez=rxz-lxz
        yez=ryz-lyz
        zez=rzz-lzz
        
        !      - nz cross derivatives 
        
        rxn=X(I,J+1,K)-X(I,J,K)
        ryn=Y(I,J+1,K)-Y(I,J,K)
        rzn=Z(I,J+1,K)-Z(I,J,K)
        rxe=((X(I+1,J+1,K)-X(I-1,J+1,K))+(X(I+1,J,K)-X(I-1,J,K)))*0.25D0
        rye=((Y(I+1,J+1,K)-Y(I-1,J+1,K))+(Y(I+1,J,K)-Y(I-1,J,K)))*0.25D0
        rze=((Z(I+1,J+1,K)-Z(I-1,J+1,K))+(Z(I+1,J,K)-Z(I-1,J,K)))*0.25D0
        
        rg13 = rye*rzn-ryn*rze
        rg23 = rxn*rze-rxe*rzn
        rg33 = rxe*ryn-rxn*rye
        
        tmp = 1.0D0/rg33
        
        rzz = DS/DSQRT((rg13*tmp)*(rg13*tmp)+(rg23*tmp)*(rg23*tmp)+1.0D0)
        rxz = (rg13/rg33)*rzz
        ryz = (rg23/rg33)*rzz
        
        lxn=X(I,J,K)-X(I,J-1,K)
        lyn=Y(I,J,K)-Y(I,J-1,K)
        lzn=Z(I,J,K)-Z(I,J-1,K)
        lxe=((X(I+1,J,K)-X(I-1,J,K))+(X(I+1,J-1,K)-X(I-1,J-1,K)))*0.25D0 
        lye=((Y(I+1,J,K)-Y(I-1,J,K))+(Y(I+1,J-1,K)-Y(I-1,J-1,K)))*0.25D0
        lze=((Z(I+1,J,K)-Z(I-1,J,K))+(Z(I+1,J-1,K)-Z(I-1,J-1,K)))*0.25D0 
        
        lg13 = lye*lzn -lyn*lze
        lg23 = lxn*lze -lxe*lzn
        lg33 = lxe*lyn -lxn*lye
        
        tmp = 1.0D0/lg33
        
        lzz = DS/DSQRT((lg13*tmp)*(lg13*tmp)+(lg23*tmp)*(lg23*tmp)+1.0D0)
        lxz = (lg13/lg33)*lzz
        lyz = (lg23/lg33)*lzz
        
        xnz=rxz-lxz
        ynz=ryz-lyz
        znz=rzz-lzz
        
        
        !en cross derivative
        xen=(X(I+1,J+1,K)-X(I-1,J+1,K)+X(I-1,J-1,K)-X(I+1,J-1,K))*0.25D0
        yen=(Y(I+1,J+1,K)-Y(I-1,J+1,K)+Y(I-1,J-1,K)-Y(I+1,J-1,K))*0.25D0
        zen=(Z(I+1,J+1,K)-Z(I-1,J+1,K)+Z(I-1,J-1,K)-Z(I+1,J-1,K))*0.25D0
        
        !P,Q,R value from RHS times Transpose of the cofactor matrix 
        !divided by Jacobian 
        
        xRHS=-JAC*JAC*(a11*xee*a22*xnn+a33*xzz+2.0D0*(a12*xen+a13*xez+a23*xnz))
        yRHS=-JAC*JAC*(a11*yee*a22*ynn+a33*yzz+2.0D0*(a12*yen+a13*yez+a23*ynz))
        zRHS=-JAC*JAC*(a11*zee*a22*znn+a33*zzz+2.0D0*(a12*zen+a13*zez+a23*znz))
        
        
        t1 = (xRHS*g11+yRHS*g21+zRHS*g31)*JAC
        t2 = (xRHS*g12+yRHS*g22+zRHS*g32)*JAC
        t3 = (xRHS*g13+yRHS*g23+zRHS*g33)*JAC
        
        !Updating the P, Q, and R values 
        
        GPJ(I,J,K)=GPJ(I,J,K)+0.03D0*(t1-GPJ(I,J,K))
        GQJ(I,J,K)=GQJ(I,J,K)+0.03D0*(t2-GQJ(I,J,K))
        GRJ(I,J,K)=GRJ(I,J,K)+0.03D0*(t3-GRJ(I,J,K))
        
        !PRINT*, 'Gp =',GPJ(I,J,K)
        !PRINT*, 'Gq =',GQJ(I,J,K) 
        !PRINT*, 'Gr =',GRJ(I,J,K) 
        !READ* 
        
        return
        end subroutine SPQR
        
        SUBROUTINE SPRDPQ(GPJ,GQJ,GRJ) 
!        include "common.h" 
        double precision :: GPJ(NIP1,NJP1,NKP1),GQJ(NIP1,NJP1,NKP1),GRJ(NIP1,NJP1,NKP1)
        INTEGER I,J,K 
        
        DO K=3,NK 
        DO J=2,NJ 
        DO I=2,NI 
         GPJ(I,J,K)=GPJ(I,J,2)*EXP(-0.5D0*(K-2)) 
         GQJ(I,J,K)=GQJ(I,J,2)*EXP(-0.5D0*(K-2)) 
         GRJ(I,J,K)=GRJ(I,J,2)*EXP(-0.5D0*(K-2)) 
        ENDDO 
        ENDDO 
        ENDDO 
        RETURN 
        END SUBROUTINE SPRDPQ
        
!	double FUNCTION ZCURV(XP,YP) 
!       double XP,YP,CURV 
!       ZCURV=
!       RETURN 
!       END FUNCTION ZCURV        
        
        SUBROUTINE UGRID(X,Y,Z) 
!        include "common.h" 
!
!THIS SUBROUTINE IS USED TO SPECIFY THE GRID TO BE USED FOR THE 
!COMPUTATION. EACH OF THE COORDINATES MUST BE SPECIFIED FOR THE 
!GRID BOUNDING THE CONTROL VOLUMES. 
        double precision :: X(NIP1,NJP1,NKP1),Y(NIP1,NJP1,NKP1),Z(NIP1,NJP1,NKP1)
        double precision ::  DELTAX,DX(NI-1),ZCURV
        double precision :: XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,XO,XP,XR,PW,RI1,RI2,RJ1,RJ2 
        double precision :: X1(NIP1,NJP1,NKP1),Y1(NIP1,NJP1,NKP1),Z1(NIP1,NJP1,NKP1) 
        double precision :: X2(NIP1,NJP1,NKP1),Y2(NIP1,NJP1,NKP1),Z2(NIP1,NJP1,NKP1)
        double precision :: DX1,TSUM,DSUM,DELTAY,Yo,Yp,DELTAZ 
        integer CN,CI,CJ 
        integer I,J,K 
        parameter (XMIN = -4.0D0, XMAX = 4.0D0) 
        parameter (YMIN = -4.0D0, YMAX = 4.0D0) 
        parameter (ZMIN = 0.0D0, ZMAX=4.0D0) 
                
        DELTAX = (XMAX-XMIN)/DBLE(NI-2)
        DELTAY = (YMAX-YMIN)/DBLE(NJ-2) 
        DELTAZ = (ZMAX-ZMIN)/DBLE(NK-2) 
        
!** Bottom boundary points redistribution to attract points to the centre of the bump using
!-- geometrical ratio XR ****
        
        CN=INT((NI-1)/2)+1
        XR=0.95D0
        
        DSUM=0.0D0 
        DO I=1,CN-2 
        CALL POW(XR,I,PW) 
        DSUM=DSUM+(1.0D0/PW) 
        PW=0.0D0 
        ENDDO
                TSUM=2.0D0*(DSUM+1.0D0) 
                DX1=(XMAX-XMIN)/TSUM 
                DX(CN)=DX1 
                DX(CN-1)=DX1 
        
        DO I=1,CN-2 
        CALL POW(XR,I,PW) 
        DX(CN+I)=DX1/PW 
        PW=0.0D0 
        ENDDO
        
        DO I=CN-2,1,-1 
        CALL POW(XR,I,PW) 
        DX((CN-1)-I)=DX1/PW 
        PW=0.0D0 
        ENDDO 
        
        K=2 
        DO I=2,NI 
        Y(I,2,K)=YMIN 
        Y(I,NJ,K)=YMAX 
        ENDDO 
        
        K=2 
        DO J=2,NJ 
        X(2,J,K)=XMIN 
        X(NI,J,K)=XMAX 
        ENDDO 
        
!--Bottom Surface 
        K=2 
        X(2,2,K)=XMIN 
        Y(2,2,K)=YMIN 
        DO I=3,NI-1 
        X(I,2,K)=X(I-1,2,K)+DX(I-2) 
        X(I,NJ,K)=X(I,2,K) 
        Y(I,2,K)=YMIN 
        Y(I,NJ,K)=YMAX 
        ENDDO 
        
        DO J=3,NJ-1 
        Y(2,J,K)=Y(2,J-1,K)+DX(J-2) 
        X(2,J,K)=XMIN 
        Y(NI,J,K)=Y(2,J,K) 
        X(NI,J,K)=XMAX 
        ENDDO 
        
        K=2 
        DO I=2,NI 
        DO J=2,NJ 
        RI1=DBLE(I-2)/DBLE(NI-2) 
        RI2=DBLE(NI-I)/DBLE(NI-2) 
        X1(I,J,K)=RI1*X(NI,J,K)+RI2*X(2,J,K) 
        Y1(I,J,K)=RI1*Y(NI,J,K)+RI2*Y(2,J,K) 
        ENDDO 
        ENDDO 
        
        DO I=2,NI 
        DO J=2,NJ 
        RJ1=DBLE(J-2)/DBLE(NJ-2) 
        RJ2=DBLE(NJ-J)/DBLE(NJ-2) 
        X2(I,J,K)=RJ1*(X(I,NJ,K)-X1(I,NJ,K))+RJ2*(X(I,2,K)-X1(I,2,K)) 
        Y2(I,J,K)=RJ1*(Y(I,NJ,K)-Y1(I,NJ,K))+RJ2*(Y(I,2,K)-Y1(I,2,K)) 
        ENDDO 
        ENDDO 
        
        DO I=2,NI 
        DO J=2,NJ 
        X(I,J,K)=X1(I,J,K)+X2(I,J,K) 
        Y(I,J,K)=Y1(I,J,K)+Y2(I,J,K) 
        ENDDO 
        ENDDO 
        
        Y(1,2,K)=Y(2,2,K)
        X(1,2,K)=XMIN-DELTAX
        Y(NI+1,2,K)=Y(NI,2,K)
        X(NI+1,2,K)=XMAX+DELTAX
        
        K=2
        DO J=2,NJ
        DO I=2,NI
        Z(I,J,K)=1.0D0*EXP(-0.5D0*(X(I,J,K)**2+Y(I,J,K)**2))  
        ENDDO 
        ENDDO 
        
!--Top Surface 
        K=NK
        XO=XMIN
        DO I=2,NI
        XP=XO+DELTAX*(I-2)
        X(I,2,K)=XP
        Y(I,2,K)=YMIN
        X(I,NJ,K)=XP
        Y(I,NJ,K)=YMAX
        ENDDO
        
        DO I=2,NI 
        YO = Y(I,2,K) 
        DO J=3,NJ-1 
        YP=YO+DELTAY 
        YO=YP 
        Y(I,J,K)=YP 
        X(I,J,K)=X(I,NJ,K) 
        X(1,J,K)=XMIN-DELTAX 
        X(NIP1,J,K)=XMAX+DELTAX 
        ENDDO 
        ENDDO 
        
        K=NK 
        DO J=2,NJ 
        DO I=2,NI 
        Z(I,J,K)=ZMAX 
        ENDDO 
        ENDDO 
        
        !--Front
        I=NI
        J=2
        DO K=2,NK
        Z(I,J,K)=ZMIN+DELTAZ*(K-2)
        X(I,J,K)=XMAX 
        Y(1,J,K)=YMIN 
        ENDDO 
        
        I=2 
        J=2 
        DO K=2,NK 
        Z(I,J,K)=ZMIN+DELTAZ*(K-2) 
        X(I,J,K)=XMIN 
        Y(I,J,K)=YMIN 
        ENDDO 
        
        J=2
        DO I=2,NI 
        DO K=2,NK 
        RI1=DBLE(I-2)/DBLE(NI-2) 
        RI2=DBLE(NI-I)/DBLE(NI-2) 
        X1(I,J,K)=RI1*(X(NI,J,K)+RI2*X(2,J,K))
        Z1(I,J,K)=RI1*(Z(NI,J,K)+RI2*Z(2,J,K))
        ENDDO 
        ENDDO 
        
        DO I=2,NI 
        DO K=2,NK 
        RJ1=DBLE(K-2)/DBLE(NK-2) 
        RJ2=DBLE(NK-K)/DBLE(NK-2) 
        X2(I,J,K)=RJ1*(X(I,J,NK)-X1(I,J,NK))+RJ2*(X(I,J,2)-X1(I,J,2)) 
        Z2(I,J,K)=RJ1*(Z(I,J,NK)-Z1(I,J,NK))+RJ2*(Z(I,J,2)-Z1(I,J,2)) 
        ENDDO 
        ENDDO 
        
        DO I=2,NI
        DO K=2,NK
        X(I,J,K)=X1(I,J,K)+X2(I,J,K) 
        Z(I,J,K)=Z1(I,J,K)+Z2(I,J,K) 
        ENDDO 
        ENDDO 
        
        J=2 
        DO K=2,NK 
        DO I=2,NI 
                Y(I,J,K)=YMIN 
        ENDDO 
        ENDDO 
        
        !--Back Surface 
        I=NI 
        J=NJ 
        DO K=2,NK 
        Z(I,J,K)=ZMIN+DELTAZ*(K-2) 
        X(I,J,K)=XMAX 
        Y(I,J,K)=YMAX 
        ENDDO 
        
        I=2 
        J=NJ 
        DO K=2,NK 
        Z(I,J,K)=ZMIN+DELTAZ*(K-2) 
        X(I,J,K)=XMIN 
        Y(I,J,K)=YMAX 
        ENDDO 
        
        J=NJ 
        DO I=2,NI 
        DO K=2,NK 
        RI1=DBLE(I-2)/DBLE(NI-2) 
        RI2=DBLE(NI-I)/DBLE(NI-2) 
        X1(I,J,K)=RI1*X(NI,J,K)+RI2*X(2,J,K) 
        Z1(I,J,K)=RI1*Z(NI,J,K)+RI2*Z(2,J,K) 
        ENDDO 
        ENDDO 
        
        DO I=2,NI 
        DO K=2,NK 
        RJ1=DBLE(K-2)/DBLE(NK-2) 
        RJ2=DBLE(NK-K)/DBLE(NK-2) 
        X2(I,J,K)=RJ1*(X(I,J,NK)-X1(I,J,NK))+RJ2*(X(I,J,2)-X1(I,J,2))
        Z2(I,J,K)=RJ1*(Z(I,J,NK)-Z1(I,J,NK))+RJ2*(Z(I,J,2)-Z1(I,J,2)) 
        ENDDO 
        ENDDO 
        
        DO I=2,NI 
        DO K=2,NK 
        X(I,J,K)=X1(I,J,K)+X2(I,J,K) 
        Z(I,J,K)=Z1(I,J,K)+Z2(I,J,K) 
        ENDDO 
        ENDDO 
        
        J=NJ
        DO K=2,NK
        DO I=2,NI
        Y(I,J,K)=YMAX 
        ENDDO 
        ENDDO 
        
        !--Left Surface 
        I=2 
        DO J=2,NJ 
        DO K=2,NK 
        RI1=DBLE(J-2)/DBLE(NJ-2) 
        RI2=DBLE(NJ-J)/DBLE(NJ-2) 
 	Y1(I,J,K)=RI1*Y(I,NJ,K)+RI2*Y(I,2,K)
       	Z1(I,J,K)=RI1*Z(I,NJ,K)+RI2*Z(I,2,K) 
        ENDDO 
        ENDDO 
        
        DO J=2,NJ 
        DO K=2,NK 
        RJ1=DBLE(K-2)/DBLE(NK-2) 
        RJ2=DBLE(NK-K)/DBLE(NK-2) 
        Y2(I,J,K)=RJ1*(Y(I,J,NK)-Y1(I,J,NK))+RJ2*(Y(I,J,2)-Y1(I,J,2))
        Z2(I,J,K)=RJ1*(Z(I,J,NK)-Z1(I,J,NK))+RJ2*(Z(I,J,2)-Z1(I,J,2)) 
        ENDDO
        ENDDO
        
        DO J=2,NJ 
        DO K=2,NK 
        Y(I,J,K)=Y1(I,J,K)+Y2(I,J,K) 
        Z(I,J,K)=Z1(I,J,K)+Z2(I,J,K) 
        ENDDO 
        ENDDO 
        
        I=2 
        DO K=2, NK 
        DO J=2, NJ 
        X(I,J,K)=XMIN 
        ENDDO 
        ENDDO 
        
        !--Right Surface 
        I=NI
        DO J=2,NJ
        DO K=2,NK
        RI1=DBLE(J-2)/DBLE(NJ-2) 
        RI2=DBLE(NJ-J)/DBLE(NJ-2) 
	Z1(I,J,K)=RI1*Z(I,NJ,K)+RI2*Z(I,2,K)
	Y1(I,J,K)=RI1*Y(I,NJ,K)+RI2*Y(I,2,K)
       	ENDDO 
        ENDDO 
        
        DO J=2 ,NJ 
        DO K=2,NK 
        RJ1=DBLE(K-2)/DBLE(NK-2) 
        RJ2=DBLE(NK-K)/DBLE(NK-2) 
        Y2(I,J,K)=RJ1*(Y(I,J,NK)-Y1(I,J,NK))+RJ2*(Y(I,J,2)-Y1(I,J,2)) 
        Z2(I,J,K)=RJ1*(Z(I,J,NK)-Z1(I,J,NK))+RJ2*(Z(I,J,2)-Z1(I,J,2)) 
        ENDDO 
        ENDDO 
        
        DO J=2,NJ 
        DO K=2,NK 
        Y(I,J,K)=Y1(I,J,K)+Y2(I,J,K) 
        Z(I,J,K)=Z1(I,J,K)+Z2(I,J,K) 
        ENDDO 
        ENDDO 
        
        I=NI 
        DO K=2,NK 
        DO J=2,NJ 
        X(I,J,K)=XMAX 
        ENDDO 
        ENDDO 
        
        !--Initial Grid
        
        DO I=3,NI-1
        DO J=3,NJ-1 
        DELTAZ=(Z(I,J,NK)-Z(I,J,2))/DBLE(NK-2)
        DO K=3,NK-1
        Z(I,J,K)=Z(I,J,K-1)+DELTAZ
        X(I,J,K)=X(I,J,2) 
        Y(I,J,K)=Y(I,J,2) 
        ENDDO 
        ENDDO 
        ENDDO 
        
        !--Save the grid to PlotMtv format file 
        CALL GRDEXT(X,Y,Z)
        CALL SAVEGRD(X,Y,Z)
        CALL SAVEMTV(X,Y,Z) 
        
        RETURN 
        end subroutine UGRID
        
        
        subroutine GRDEXT(X,Y,Z) 
!        include "common.h" 
        double precision :: X(NIP1,NJP1,NKP1) 
        double precision ::  Y(NIP1,NJP1,NKP1) 
        double precision :: Z(NIP1,NJP1,NKP1) 
        double precision :: tmp1 
        integer I,J,K 
        
        DO K=2,NK
        DO J=2,NJ
        X(1,J,K)=2.0D0*X(2,J,K)-X(3,J,K) 
        X(NIP1,J,K)=2.0D0*X(NI,J,K)-X(NI-1,J,K)
        Y(1,J,K)=Y(3,J,K)
        Y(NIP1,J,K)=Y(NI-1,J,K) 
        Z(1,J,K)=Z(3,J,K) 
        Z(NIP1,J,K)=Z(NI-1,J,K) 
        ENDDO 
        ENDDO 
        
        DO K=2,NK 
        DO I=2,NI 
        X(I,1,K)=X(I,3,K) 
        X(I,NIP1,K)=X(I,NJ-1,K) 
        Y(I,1,K)=2.0D0*Y(I,2,K)-Y(I,3,K) 
        Y(I,NIP1,K)=2.0D0*Y(I,NJ,K)-Y(I,NJ-1,K) 
        Z(I,1,K)=Z(I,3,K) 
        Z(I,NJP1,K)=Z(I,NJ-1,K) 
        ENDDO 
        ENDDO 
        
        !---Corner Points 
        
        DO K=2,NK 
        tmp1=2.0D0*X(2,1,K)-X(3,1,K) 
        X(1,1,K)=2.0D0*X(2,2,K)-X(3,3,K) 
        tmp1=2.0D0*X(2,NJP1,K)-X(3,NJP1,K) 
        X(1,NJP1,K)=2.0D0*X(2,NJ,K)-X(3,NJ-1,K) 
        tmp1=2.0D0*X(NI,1,K)-X(NI-1,1,K) 
        X(NIP1,1,K)=2.0D0*X(NI,2,K)-X(NI-1,3,K) 
        tmp1 = 2.0D0*X(NI,NJP1,K)-X(NI-1,NJP1,K) 
        X(NIP1,NJP1,K)=2.0D0*X(NI,NJ,K)-X(NI-1,NJ-1,K) 
        
        tmp1 = 2.0D0*Y(1,2,K)-Y(1,3,K) 
        Y(1,1,K)= 0.5D0*(tmp1+0.5D0*(Y(2,1,K)+Y(3,1,K))) 
        tmp1 = 2.0D0*Y(1,NI,K)-Y(1,NI-1,K) 
        Y(1,NJP1,K)= 0.5D0*(tmp1+0.5D0*(Y(2,NJP1,K)+Y(3,NJP1,K))) 
        tmp1 = 2.0D0*Y(NIP1,2,K)-Y(NIP1,3,K) 
        Y(NIP1,1,K)= 0.5D0*(tmp1+0.5D0*(Y(NI,1,K)+Y(NI-1,1,K))) 
        tmp1 = 2.0D0*Y(NIP1,NI,K)-Y (NIP1,NI-1,K) 
        Y(NIP1,NJP1,K)=0.5D0*(tmp1+0.5D0*(Y(NI,NJP1,K)+Y(NI-1,NJP1,K))) 
        
        Z(1,1,K)=Z(3,3,K) 
        Z(1,NJP1,K)=Z(3,NJ-1,K) 
        Z(NIP1,1,K)=Z(NI-1,3,K) 
        Z(NIP1,NJP1,K)=Z(NI-1,NJ-1,K) 
        enddo
        return 
        end subroutine GRDEXT
        
        SUBROUTINE SEESPC(X,Y,Z)
!        include "common.h"
        double precision :: X(NIP1,NJP1,NKP1),Y(NIP1,NJP1,NKP1),Z(NIP1,NJP1,NKP1)
        double precision :: ADS,DX,DY,DZ,MAXADS
        integer I,J
        MAXADS=0.0D0
        DO J=2,NJ
        DO I=2,NI
        DX=X(I,J,3)-X(I,J,2)
        DY=Y(I,J,3)-Y(I,J,2) 
        DZ=Z(I,J,3)-Z(I,J,2) 
        ADS=DSQRT(DX*DX+DY*DY+DZ*DZ) 
        MAXADS=MAX(MAXADS,ADS) 
        ENDDO 
        ENDDO 
        PRINT*,'MAX DS  ',MAXADS 
        return
        end subroutine SEESPC
        
        SUBROUTINE POW(X,N,PW) 
        double precision :: X,TMP,PW 
        INTEGER N 
        TMP=N*LOG(X) 
        PW=EXP(TMP) 
        RETURN 
        END SUBROUTINE
        
        SUBROUTINE SAVEMTV(X,Y, Z) 
!        include "common.h" 
        double precision :: X(NIP1,NIP1,NKP1) 
        double precision :: Y(NIP1,NIP1,NKP1) 
        double precision :: Z(NIP1,NIP1,NKP1) 
        INTEGER I,J,K 
        
        OPEN(23,FILE='3dgrid.mtv') 
        WRITE(23,998)'$ DATA = CURVE3D' 
        WRITE(23,998)'% axisscale = False' 
        WRITE(23,998)'% axislabel = False' 
        
        I=16
        DO K=2,NK
        DO J=2,NJ
        WRITE(23,*) X(I,J,K),' ',Y(I,J,K),' ',Z(I,J,K)
        ENDDO 
        WRITE(23,*) 
        ENDDO 
        WRITE(23,*) 
        
        DO J=2 ,NJ 
        DO K=2,NK 
        WRITE(23,*) X(I,J,K),' ',Y(I,J,K),' ',Z(I,J,K) 
        ENDDO 
        WRITE(23,*) 
        ENDDO 
        WRITE(23,*) 
        
        J=16 
        DO K=2,NK 
        DO I=2,NI 
        WRITE(23,*) X(I,J,K)," ",Y(I,J,K)," ",Z(I,J,K)
        ENDDO 
        WRITE(23,*) 
        ENDDO 
        WRITE(23,*) 
        
        K=2
        DO J=12,NJ
        DO I=16,NI
        WRITE(23,*) X(I,J,K),' ',Y(I,J,K),' ',Z(I,J,K)
        ENDDO
        WRITE(23,*) 
        ENDDO 
        WRITE(23,*) 
        
        DO I=16,NI
        DO J=2,NJ
        WRITE(23,*) X(I,J,K),' ',Y(I,J,K),' ',Z(I,J,K) 
        ENDDO 
        WRITE(23,*) 
        ENDDO 
        WRITE(23,*) 
        
        CLOSE (23) 
998     FORMAT(A) 
        return
        end subroutine SAVEMTV
        
        SUBROUTINE SAVEGRD(X,Y,Z) 
!        include "common.h" 
        double precision :: X(NIP1,NJP1,NKP1),Y(NIP1,NIP1,NKP1),Z(NIP1,NIP1,NKP1) 
        INTEGER I,J,K 
        OPEN(24,FILE='3dgrid.plt') 
        WRITE(24,*)'VARIABLES = "X","Y","Z","T"'
        WRITE(24,*)'ZONE F=POINT,I=',NI-1,',J=',NJ-1,',K=',NK-1 
        DO K=2,NK 
        DO J=2,NJ 
        DO I=2,NI 
        WRITE(24,*) X(I,J,K),' ',Y(I,J,K),' ',Z(I,J,K),'0.0'
        ENDDO 
        ENDDO 
        ENDDO 
        WRITE(24,*) 
        CLOSE(24) 
1001    FORMAT(A) 
        return 
        end subroutine SAVEGRD        
        
        
!       END OF MAIN PROGRAM        
	END PROGRAM XSTREAM
